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ABSTRACT 


A  simple  model  of  a  time  sequence  of  star  images  containing  a  moving  point  object 
(satellite)  is  developed.  Optimal  signal  enhancement  and  detection  processing  theory  is  iq>- 
plied  to  this  model  and  a  three-dimensional  Fourier  matched  filter  implementation  is  derived 
to  compute  clutter-to-noise  ratio  (CNR)  suppression,  signal-to-noise  ratio  (SNR)  enhance¬ 
ment,  and  probability  of  detection  (PJ)  and  false  alarm  (Pj^)  rate  estimates  as  a  function  of 
input  single  pixel  SNR.  Using  this  theory  allows  one  to  compute  the  Fourier  domain  matched 
filters  directly,  thereby  eliminating  the  enormous  storage  cost  associated  with  large  banks  of 
three-dimensional  matched  filters.  The  model  and  theory  are  tested  using  computer-generated 
simulated  data  sets  having  known  noise  and  clutter  characteristics.  The  theory  is  then  applied 
to  real  data  sets  collected  using  a  Lincoln  Laboratory  CCD  camera.  Initial  results  indicate 
good  agreement  with  matched  filter  theory.  Detection  of  objects  with  an  initial  single  pixel 
SNR  B  1  is  demonstrated. 
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1.  INTRODUCTION 


Satellite  detection  from  ground-  and  space-based  optical  sensors  entails  acquiring  a  sequence  of 
two-dimensional  intensity  images,  thereby  forming  a  three-dimensional  image  volume  with  two  spatial 
coordinates  and  one  time  coordinate.  Lincoln  Laboratory  uses  a  camera  containing  a  CCD  chip  with  420 
X  420  pixels  in  the  focal  plane.  Versions  of  the  chip  are  available  that  can  be  mounted  flush  with  each 
other  to  give  840  X  840  or  mote  pixels  per  frame.  Most  of  the  existing  data  sets  used  in  this  study  consist 
of  16  consecutive  frames  of  data,  although  many  more  could  easily  be  collected.  Integration  times  per 
frame  ate  typically  0.2  to  5  s.  The  result  of  this  data  collection  process  is  an  image  volume  of  x 
X  N  pixels,  where  N  and  N  are  the  number  of  pixels  in  the  two-dimensional  CCD  focal  plane  array, 
and  N,  is  the  number  of  frames  of  data. 

Images  collected  in  this  process  typically  contain  noise,  stars,  and  possibly  one  or  more  satellites. 
The  stars  and  the  satellites  are  imaged  as  unresolved  point  sources.  The  most  common  method  of  data 
collection  is  sidereal  track  mode,  wherein  the  telescope  tracks  the  stars,  which  remain  fixed  on  the  focal 
plane.  In  the  sidereal  case  a  satellite  moves  across  the  detector,  forming  a  steak.  Another  data  collection 
method  is  the  ephemeris  track  mode,  which  requires  an  estimate  of  the  satellite  orbit.  With  this  method 
the  telescope  tracks  the  satellite,  which  remains  fixed  on  the  focal  plane,  while  the  stars  form  streaks  on 
the  detector.  To  detect  a  satellite  a  differentiation  must  be  made  between  the  stars  (moving  at  one  velocity 
on  the  focal  plane)  and  a  satellite  (moving  at  a  different  velocity).  Furthermore,  in  the  case  of  a  faint 
satellite  the  single  pixel  signal-to-noise  ratio  (SNR)  may  be  considerably  less  than  one,  making  it  essen¬ 
tially  undetectable  to  the  unaided  observer.  Thus  there  are  two  tasks  to  be  performed:  clutter  (star) 
suppression  and  signal  enhancement,  shown  schematically  in  Figure  1. 
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Figure  I.  Matched  filter  detection  processing  chain. 
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The  purpose  of  this  report  is  to  develop  and  test  a  model  of  the  image  sets  used  for  satellite 
detection  and  a  corresponding  detection  theory.  Insights  gained  from  this  work  are  used  to  evaluate 
current  state-of-the-art  signal  processing  algorithms.  Probability  of  detection  (Pj)  and  of  false  alarm 
curves  are  determined  for  the  theory  as  a  function  of  single  pixel  SNR.  Theoretically  predicted  SI^ 
enhancement  values  can  then  be  compared  with  results  from  real  data  sets  to  assess  the  validity  of  the 
current  model.  The  discussion  begins  with  matched  filter  theory  as  the  optimal  linear  technique  to 
determine  the  presence  and  location  of  a  known  signal  in  a  background  of  stationary  noise  [1].  If  the  input 
signal  is  sufficiently  well  behaved  (stationary  Gaussian  statistics),  the  linear  matched  filter  represents  the 
optimal  processing  at  any  order.  Detection  can  be  performed  using  either  a  correlation  function  in  the 
space-time  volume  or  the  mathematically  equivalent  matched  filter  approach  in  the  Fourier  transform 
volume  [2].  Reed  [3]  has  extended  the  common  one-dimensional  filter  used  in  radar  applications  and  the 
two-dimensional  filter  used  for  images  to  a  three-dimensional  time  sequence  of  images.  Reed  [4]  and  Porat 
[S]  have  each  discussed  the  application  of  matched  filters  to  the  three-dimensional  point  object  detection 
problem.  The  nature  of  the  Fourier  volume  for  this  application  and  the  form  of  the  matched  filter  will 
be  explored.  Next  is  the  question  of  detectability:  What  is  the  relation  between  satellite  brightness,  noise 
level,  observation  time,  and  detection  capability? 

The  theory  is  initially  developed  for  a  continuous  case  (both  spatial  and  temporal)  for  analytic 
simplicity  and  clarity.  (A  discrete  formulation  is  presented  in  the  appendix,  and  data  representing  a 
discretely  sampled  image  volume  is  generated  and  used  to  verify  the  theory.)  To  perform  the  signal 
enhancement  step,  one  must  first  perform  clutter  suppression.  Thus  three  methods  are  explored:  a  fre¬ 
quency  slice  based  on  the  matched  filter  theory,  a  maximum  likelihood  normalizer,  and  a  star  map 
technique  that  introduces  additional  information  to  suppress  the  clutter.  Results  for  two-  and  three- 
dimensional  data  sets,  both  computer  generated  and  real,  are  presented.  These  data  sets  verify  the  algorithm 
and  allow  a  comparison  of  the  model  with  real-life  data,  which  introduce  Poisson  noise,  mount  jitter,  and 
effects  due  to  the  finite  size  of  the  CCD  cells. 
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2.  THEORY 


2.1  THE  SPATIAL  DOMAIN  MODEL 

To  understand  the  image  data  and  generate  a  sensitive  detection  aigorithm,  it  is  useful  to  form  an 
analytical  model  of  the  image  data.  Scribner  has  developed  a  model  for  the  entire  image  chain,  from 
object  radiance  through  the  imaging  hardware  to  the  final  electronic  image  [6].  This  model  includes 
factors  for  the  optics  point  spread  function  (PSF)  and  the  effects  of  discrete  sampling.  A  simpler  approach 
is  taken  here  with  the  focus  on  the  image  set.  The  signal  intensity  at  any  given  pixel  will  be  the  sum  of 
a  white  noise  contribution  (7^  and  possible  contributions  from  stars  (/^)  and/or  satellites  (/p,  that  is,  I(jc,/)  = 
/^(ac,f)  +  I^x,t)  + 

The  discrete  image  volume  is  composed  of  N^y.  pixels,  each  having  volume  =  Pj^PyP,- 
The  total  image  volume  is  then  V,  =  LJLyT,  where  =  P^x'  ~  P/^y'  ^  ~  P^t  convenience 

is  denoted  by  (L,T). 

Both  the  stars  and  satellites  may  be  described  as  diffraction-limited  Gaussian  spots  falling  on  the 
focal  plane  of  a  telescope  with  the  same  Gaussian  width  parameter  o.  For  now  the  Poisson  nature  of  the 
photon  statistics  is  ignored  and  constant  irradiance  is  assumed.  Stars  are  described  by  the  intensity 
distribution 


(I) 


where  the  sum  is  over  all  stars  and  they  are  moving  at  the  same  velocity  v.:x.(r)  =  x.  +  v.r.  A  satellite 
in  the  image  is  described  by  a  similar  distribution. 


I^{x,t)  =  A^it)e  .  (2) 

The  satellite  is  moving  at  a  different  velocity,  vy.x^(i) v/-  The  velocities  are  two-dimensional  focal 
plane  velocities  (e.g.,  millimeters  per  second).  In  discrete  focal  plane  coordinates  could  be  written  (in 
one  dimension)  as  =  Px^^Pv  pixels/frame.  The  time  dependence  of  AJj)  represents 

tumbling  or  rotating  objects  as  well  as  changes  in  the  sun,  satellite,  or  observer  angle.  For  this  study  the 
time  independent  form  Ap)  =  A^  is  used.  The  units  of  A^  and  the  A  -  are  irradiance: 


^  r  j  energy  ..power 
^  area  •  time  ^  area 


The  rest  of  this  section  discusses  only  the  signal  due  to  a  satellite;  howevi  '  scause  the  PSFs  are 
identical,  the  stars  behave  identically. 
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The  energy  deposited  per  frame  time  (p,)  by  a  satellite  is 


Pi  +CO 

0  —  OO 

0  — ®o 

EXp,)=2„a^A^p,  .  (3) 

While  the  focal  plane  is  obviously  finite,  it  is  sufficiently  large  compared  with  the  Gaussian  PSF  to  use 
the  infinite  integrals  for  all  streaks  except  those  that  actually  leave  the  focal  plane.  This  energy  may  be 
confined  to  a  single  pixel  or  spread  among  several,  depending  on  v^. 

A  useful  reality  check  after  the  transforms  are  made  is  Parseval’s  Theorem  P  =  lf\x)dx  =  ig^(ai)d(o: 

P  =  J  J  . 

J  J  J  ^  ^  area*time 

0  — 

An  average  energy  deposited  per  area  along  the  streak  in  the  focal  plane  is  easily  computed  by 

g  -  energy  /  frame 
^  streak  length* streak  width 

If  the  Gaussian  spot  size  is  defined  as  2a,  the  streak  length  may  be  written  as  2a+l  v^lp^,  giving 


e 

s 


e 

s 


£ 

U'’.) 

|2a+ 

V 

s 

1 

Pt. 

•  2a 

2a+ 


K\P, 


However,  to  compute  a  single  pixel  SNR  an  energy  deposited  per  pixel  must  be  computed.  In  general, 
even  with  an  infinitesimal  point  distribution,  a  single  value  is  not  observed  because  a  point  moving  at 
an  angle  with  respect  to  the  two  spatial  axes  moves  diagonally  through  the  pixel  grid,  thereby  cutting 
through  comers  of  some  pixels  and  depositing  more  energy  in  some  than  in  others.  The  situation  is  even 
more  complicated  when  using  a  finite  PSF. 

If  the  physical  dimensions  of  the  CCD  pixels  are  and  the  spot  area  (27ia^)  is  much  less  than 

the  pixel  area,  then  for  a  streak  parallel  to  the  jr-axis  and  centered  on  the  row  of  pixels,  two  general  cases 
can  be  distinguished.  If  the  velocity  is  sufficiently  slow  that  the  object  can  ly*  entirely  contained  within 
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a  single  pixel,  then  the  energy  per  pixel  is  the  energy  per  frame.  Note  however,  that  as  long  as  v  ^  0, 
there  are  some  frames  where  the  energy  is  distributed  between  two  pixels.  If  the  velocity  times  integration 
time  is  greater  than  a  single  pixel,  then  the  energy  per  pixel  is  the  energy  per  frame  times  the  ratio  of 
a  pixel  length  to  the  streak  length.  Thus  the  energy  per  pixel  is 


e  = 

S 


(4) 


which  will  be  used  as  a  general  approximation  to  the  single  pixel  SNR.  It  is  also  useful  to  determine  the 
total  number  of  pixels  contained  in  such  a  streak: 


Equation  (S)  undercounts  the  true  number  of  pixels  in  a  streak  because  of  the  discrete  nature  of  the 
frames.  A  better  formula  would  compute  I  /p^,  round  up  to  the  next  highest  integer,  and  multiply 
by  N^.  For  this  work  the  simpler  formulation  is  sufficient 

Noise  in  the  image  volume  may  be  described  by  a  white  (in  both  space  and  time)  Poisson  dailc 
current,  which  may  be  written  as  a  constant  component  Jf  plus  a  zero  mean  noise  component  of  amplitude 
£:Ti=ri±£,  where  e  is  drawn  from  a  normal  distribution  characterized  by  standard  deviation  p  =  yfrj. 
For  the  fluctuation  portion  the  autocorrelation  function  is  written 

,  (6) 

and  the  expectation  value  in  any  volume  element  (pixel)  is  written 

(/,(i.0)=5v,  =  l, 

where  pp  is  the  observed  standard  deviation  of  the  pixels.  The  total  noise  energy  is  then  ,but  the 

Parseval  “energy”  is  P=  +p^V  =  (rj^  +  Tf)V  . 

Other  phenomena  that  are  observed  in  actual  images  but  ignored  in  this  model  include 

•  Fixed  pattern  noise 

•  Saturation  effects 

•  Time  varying  effects  resulting  from  cloud  motion 

•  Atmospheric  and  mount  jitter 
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•  Radiation  effects  (gamma  events) 

•  Target  signature 

•  Photon  Poisson  statistics  of  stars  and  targets. 

Z2  DETECTION  THEORY 

Several  optimal  processing  algorithms  are  available,  depending  on  whether  the  objective  is  to 

•  Determine  if  a  signal  of  known  signature  is  present  in  a  random  noise  background,  and  if  it  is, 
to  locate  it  (linear  matched  filter) 

•  Determine  if  a  random  signal,  having  known  power  spectra,  is  present  in  a  random  noise  back¬ 
ground  (stochastic  matched  filter) 

•  Recover,  as  well  as  possible,  the  signature  of  an  unknown  signal  in  a  random  noise  background 
(Wiener  or  least  mean  square  filter). 

Each  task  is  slightly  different  and  requires  a  different  mathematical  approach  [7,S]. 

For  a  simple,  one-dimensional,  time  varying  signal  s(r),  its  Fourier  transform  (amplitude  spectrum) 
can  be  written  as  S((o)  =  An  additive  noise  source  r](/)  with  autocorrelation  function 

Rft)  =  ^n(^)i1(^-'C))  has  a  Fourier  power  spectrum  IN(q>)I^  =  !F(R(t)).  In  this  case  the  optimal  detection 
matched  filter  is  [1] 


//((»)  = 


S*(0})e  " 
IN(o))|2 


where  5* (to)  is  the  complex  conjugate  of  S(<o).  The  phase  factor  allows  the  output  detection  peak  to  be 
shifted  from  the  true  object  location  x  to  x  +  r^.  In  practice  r^  is  usually  set  to  zero.  The  filtering  operation 
is  described  s'(r)  =  !r\S(<o)ff(<o)).  Detections  now  correspond  to  peaks  in  the  output  signal  s'(t).  A 
threshold  can  be  set  at  a  given  level  above  the  background  noise  level  to  yield  a  desired  Pj  or  Pj^.  If  the 
signal  is  a  rectangle  function,  then  filtering  with  the  above  matched  filter  gives  a  triangle  function  output: 
an  optimal  detection  but  bearing  little  resemblance  to  the  input  signal.  In  three  dimensions  this  detection 
becomes 


WA 

Because  the  matched  filter  is  designed  to  optimize  the  output  SNR,  the  well-known  result  of  using 
a  matched  filter  in  the  presence  of  white  noise  (N(k,(0)  =  N)  is  that  the  final  SNR  depends  only  on  the 
total  signal  energy  and  not  its  spatial  distribution.  For  the  current  situation,  when  the  signal  increases 
linearly  with  observation  time  and  the  noise  increases  as  the  square  root,  the  matched  filter  is  expected 


(7) 
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to  increase  the  SNR  by  the  square  root  of  the  number  of  frames  of  data  collected.  The  signal  enhancement 
process  is  depicted  geometrically  in  Figure  2. 

Referring  to  Figure  2(a),  for  the  objects  that  are  present  for  the  entire  observation  time,  the  signal 
peaks  in  a  predetermined  plane.  In  this  case  one  need  only  back  transform  the  single  desired  temporal 
plane  to  perform  the  detection.  This  procedure  is  accomplished  following  the  filtering  operation  by 
performing  the  inverse  temporal  transform  only  for  t  =  however,  if  one  desires  to  perform  further 
operations  in  addition  to  basic  detection,  such  as  position  estimation,  the  use  of  the  full  three-dimensional 
result  may  be  required. 

The  result  of  performing  a  two-dimensional  back  transform  is  a  single  x-y  frame  of  data  having  a 
correlation  spike  at  the  r  =  location  of  the  target.  The  shape  of  this  spike  is  the  autocorrelation  of  a 
rectangle  the  length  of  a  single  frame  streak,  whereas  the  full  three-dimensional  form  produces  the 
autocorrelation  of  the  entire  streak.  Two  things  are  accomplished:  first,  the  size  of  the  back  transform  that 
must  be  applied  for  each  filter  is  reduced,  and  second,  by  coll^sing  the  broad  triangle  to  the  much 
narrower  single  frame  correlation  triangle,  detection  processing  is  simplified. 

The  matched  filter  operation  in  the  Fourier  volume  corresponds  to  the  correlation  function  in  the 
space-time  domain: 


+00 

ir-'(F(fi))G*(fl)))=  J/(T  +  /)g(T)dT 


In  fact,  for  the  filter  H{(a)  =  5* (to),  the  matched  filter  is  an  autocorrelation  function.  The  Fourier  domain 
arbitrary  phase  factor  is  equivalent  to  the  addition  of  a  time  delay  in  the  argument  to  the  filter  function. 

In  three  dimensions  the  enhanced  signal  is 


=  8iXf'^Mx  +  x,r+t-t^)dxdr 

+00 

To  illustrate,  using  a  5-function  target  moving  at  constant  velocity  v  observed  for  the  interval 
0  <  T<  T: 

,  T)  = /(  t)5( J  -  VT) 
and  use  the  following  relation. 


J  S(x)S(x-i-x)dX  =  S(x) 
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Figure  2.  Signal  enhancement  techniques:  (a)  Fourier  transform  matched  filter  processing  and 
(b)  space-time  convolution  processing,  "shift-and-add. " 


the  correlation  reduces  to  a  one-dimensional  integral: 


+«• 

f(x,t)=  J  Iir)I(T+t-t^)S(x-x^-vT)dT 


(8) 


The  5  function  serves  to  select  the  proper  jc  position  as  a  function  of  time.  Because  a  priori  knowledge 
of  the  target  signature  /(T)  is  lacking,  let  the  filter  profile  become  a  rectangle  function  of  amplitude  /^. 
If  the  spatial  dependence  is  temporarily  dropped  for  simplicity.  Equation  (8)  becomes: 


/(0  = 


t+T 


T+l 

o 

j  I^I(r)dx 


t 

0 


-7  +  f  <t<t 
o  o 


t  <t<T+t 

o  o 

otherwise 


(9) 


This  equation  is  a  useful  form  for  exploring  the  nature  of  the  enhanced  signal.  In  the  case  of  a 
constant  amplitude  signal  (/(r)  =  /^),  the  result  is  a  triangle  of  width  27  and  amplitude  7^7  centered  at 


/(r)  = 


/2(f+7-r  ) 

o' 

/2(r-^7-^  ) 

o'  o' 


0 


-7-t-f  <t<t 
o  o 

t  <t<T+t 
o  o 

otherwise 


If  I(t)  remains  strictly  positive,  the  result  will  still  be  centered  at  but  no  longer  symmetric.  For 
a  constant  brighmess  target  in  the  presence  of  noise,  /(/)  can  be  replaced  with  7^  -i-  h{t),  where  h{t)  is 
a  normally  distributed  random  variable  characterized  by  width  r.  The  initial  single  point  SNR  is  then 
simply  IJhit)  or  on  average  I  Jr.  The  enhanced  signal  at  /  =  is  IJT  +  «  7^7  +  7 which 

represents  an  SNR  increase  of  Vt. 
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23  FOURIER  VOLUME  DEVELOPMENT 


The  extent  of  the  Fourier  volume  is  determined  by  the  resolution  in  the  image  volume  while  the 
resolution  in  the  Fourier  volume  is  determined  by  the  extent  of  the  image  volume: 

L  LT  V 
1/-'  rr  y  -  ^ 

^  {iJtf  {2jtf 

P.PyP. 

^  (2;r)3  (2nf 

For  analytical  purposes  it  is  convenient  to  change  the  origin  to  7/2  so  that  integration  is  from  -7/2  to 
-1-7/2  rather  than  0  to  7.  The  three-dimensional  Fourier  transform  is: 


G(i,a))  =  7(/(x,/))  =  r^l  I  dt 

'-■y  -7/2  — 

-,3+7'/2  +  «r  ,  __ 


=  Gjj  (* ,  m)  -f  (* ,  fli)  -i-  (* ,  m) 

Using  the  current  model  for  the  detected  image  volume,  the  Fourier  amplitude  spectrum  may  now 
be  computed.  The  Fourier  spectrum  for  a  satellite  is 


-►r/2+«,  (-^-y*)) 

=  ^  HJV  e-^-^-‘-'dxd, 


-7/2 -< 


- -  /  '  ' 

-7/2 


_  _  a^k»k 
a^A  7e"***^e  2 


■sinc[(o>  +  **v,)j] 


(10) 
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Lacking  knowledge  of  the  expected  target  location,  the  matched  filter  will  be  the  complex  conjugate  of 
the  position-independent  portion  of  Equation  (10).  which  is  strictly  real.  Hence,  //(to)  =  ff*(0)). 


a^k»k 


- sinc|^(fl)  +  * 


(11) 


An  equivalent  representation  using  the  discrete  Fourier  transform  is  presented  in  the  Appendix.  At  this 
point  it  is  worth  noting  that  the  streak  in  the  space-time  domain  represented  the  convolution  of  a  Gaussian 
PSF  with  an  infinitesimal  line.  Thus  as  expected,  the  previous  expression  is  the  Fourier  transform  of  the 
Gaussian  multiplied  by  the  Fourier  transform  of  the  streak,  according  to  the  identity 

=  (12) 

The  Fourier  power  spectrum  is 


S _ 

2n 


sinc^ 


(13) 


If  o  0  such  that  a^A^  -»  B^,  then 


The  signal  is  maximized  in  the  k 


=  -to  plane: 


rfl27’2  g^j- 

=  — - = 

2k  So> 


B^T 


where  the  smallest  resolvable  gradation  in  the  frequency  domain,  5(o  =  2ji/T,  has  been  used  and  is  called 
the  width  of  the  plane.  As  observation  time  increases  more  signal  is  collected,  while  at  the  same  time 
occupying  a  smaller  volume  in  Fourier  space.  The  preceding  discussion  demonstrates  that  the  contribution 
due  to  the  satellite  is  confined  to  a  single  plane,  k  •  v^=  -to,  and  as  o  ->  0  it  is  uniformly  distributed  in 
that  plane  (no  k,0D  dependence). 

As  a  check  on  these  computations,  Parseval’s  theorem  is  invoked: 


2k 


-sine* 


dkdm 
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For  simplicity,  set  v  =  0  so  that 
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The  noise  Fourier  power  spectrum  will  have  a  constant  component  due  to  the  Poisson  fluctuations, 
as  well  as  a  zero  frequency  spike  corresponding  to  the  space-time  average  value,  T|.  The  fluctuations 
contribute 


Integrating  over  the  Fourier  volume  gives  P-rj^V^+p^V^  as  in  the  space-time  domain.  Characteristic 
of  white  noise  is  that  the  power  spectrum  is  a  constant  (except  for  the  zero  frequency  component) 
independent  of  location  in  the  Fourier  volume.  Also  noted  is  that  while  the  Fourier  signal  power  com¬ 
ponents  rise  as  7^,  the  average  noise  contribution  only  rises  as  T. 


2.4  CLUTTER  SUPPRESSION 

Before  performing  the  signal  enhancement  operation,  one  must  eliminate  background  clutter.  The 
optimal  method  is  a  maxii .  n  likelihood-derived  normalizer.  Because  a  detailed  discussion  and  results 
are  given  elsewhere  by  Pohlig  [9],  this  report  will  only  present  some  brief  results  using  the  normalizer 
as  a  front  end  clutter  suppression  algorithm.  The  normalizer  determines  the  average  and  standard  devia¬ 
tion  of  the  pixels  at  each  spatial  location  (a  pixel-by-pixel  average  over  time)  and  then  for  each  frame 
replaces  the  value  by  [value  -  average/standard  deviation]^. 

I.  ,-l  . 

j  _  i,j.k  t<J 

'.7.*  a 

‘•7 

This  procedure  removes  any  fixed  patterns  in  the  image  (stars,  fixed  pattern  noise,  etc.).  Because  the  stars 
represent  statistically  nonstationary  noise  (due  to  their  larger  Poisson  contribution),  the  normalizer  has  the 
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desirable  effect  of  restoring  stationary  statistics.  However,  if  one  uses  all  N,  time  samples  to  define  the 
mean  and  standard  deviation,  the  resulting  data  exhibit  a  ^  distribution  [10],  while  if  one  excludes  the 
current  temporal  value  and  uses  the  remaining  N,-l  samples,  the  data  exhibit  an  F  distribution.  If  one  uses 
the  F  distribution  approach,  the  “normalized”  data  are  in  fact  simply  a  map  of  the  single  pixel  SNR. 


Alternatively,  stars  can  be  thought  of  as  a  structured  noise  component  to  be ‘removed  by  the 
matched  filter.  In  this  case,  the  IN((a)P  term  in  the  matched  filter  description  contains  both  the  white 
background  noise  and  the  star  contribution: 


f 


Using  this  description  in  the  matched  filter  essentially  notches  out  the  intersection  of  the  spectral  plane 
of  the  stars  with  that  of  the  satellites.  The  equations  of  the  two  planes  are 


k  •v^  =-0) 

k  •V  =-Q)  , 

j 

and  the  line  of  intersection  of  the  two  planes  computed  is 


k 


y 


k 


X 


These  computations  can  efficiently  separate  a  moving  object  from  stars  that  are  moving  at  a 
different  velocity.  (This  technique  works  equally  well  in  ephemeris  track  or  sidereal  mode.)  In  practice, 
one  would  simply  remove  the  star  plane  and  then  test  with  the  candidate  velocity  filters  rather  than 
weighting  each  velocity  filter  for  the  intersection  line.  Unfortunately,  a  bright  star  located  on  a  pixel 
boundary  and  subject  to  small  position  fluctuations  generates  large  random  intensity  variations  that  are 
spread  throughout  the  Fourier  volume  and  consequently  cannot  be  removed.  Even  if  no  jitter  were 
present,  the  stars  would  still  exhibit  random  intensity  fluctuations  due  to  the  Poisson  nature  of  the  photon 
arrival  statistics. 

The  use  of  a  star  catalog  combined  with  the  telescope  line  of  sight  to  produce  a  star  mask  can 
remove  stars  to  the  14th  magnitude.  This  approach  is  extremely  powerful  and  overcomes  the  limitations 
of  the  preceding  method  by  introducing  additional  information.  Because  one  sets  the  area  of  the  star  to 
the  average  background,  minor  star  motion  and  resulting  fluctuations,  as  well  as  Poisson  noise,  can  be 
completely  eliminated  without  significantly  disturbing  the  noise  statistics.  Unlike  the  previous  method 
however,  the  star  map  does  not  remove  any  fixed  pattern  noise.  If  this  noise  presented  a  problem,  the  star 
map  could  easily  be  combined  with  a  zero  frequency  notch  filter  to  add  fixed  pattern  suppression. 
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2^  ANALYSIS  AND  DETECTION 


The  dark  current  white  noise  is  uniformly  distributed  throughout  the  Fourier  volume  (except  for  the 
zero  frequency  spike).  The  signal  due  to  stars  is  maximized  in  the  4  •  -<o  plane,  while  the  signal  from 

a  satellite  is  maximized  in  the  k  •  =  -<o  plane.  In  both  cases  the  signal  is  maximized  in  the  indicated 

plane  and  falls  off  as  a  sine  function  on  both  sides.  The  longer  the  data  collection  time,  the  better  the 
signal  is  confined  to  the  plane.  For  an  infinite  streak  the  signal  is  contained  entirely  within  a  single  plane, 
while  for  a  short  streak  the  sine  function  can  have  significant  power  in  some  of  the  lobes.  An  infinitesimal 
cross  section  for  a  streak  (single  pixel  width)  leads  to  uniform  intensity  in  the  main  signal  plane.  As  the 
signal  smears  out,  either  due  to  saturation  of  the  detector,  atmospheric  jitter,  or  resolution  of  an  object, 
the  intensity  distribution  in  the  plane  loses  its  uniform  characteristic. 

While  it  is  fairly  obvious  that  the  signal  is  maximized  in  the  k»v^=-ai  plane,  there  are  sidelobes; 
it  would  be  useful  to  know  their  amplitude  and  extent.  In  general,  the  signal  in  the  Fourier  volume  is 


c,(i-,«)|  - 


sinc^  Jfc  •v  +  0)) 


f]  ■ 


As  long  as  a\k\«  I,  the  functional  dependence  is  entirely  within  the  sine  function.  For  a  discretely 
sampled  image  volume  let  it  ->  mlnIL,  <o  ->  rjln/T,  thus  rewriting  the  sine  function  (and  defining  Sk  =  2ji/L, 
6o)  =  2n/T): 


.  ( m^vT  nScoT 

sin  - + - 

\  2  2 

mSkvT  nScoT 

- + - 

2  2 


')  .  f  mnvT  ')  .  ( m7CvT\ 

j  sin|^  ^  ng  j  sin^  ^  j 


mrevT 


mitvT 


which  is  not  very  interesting;  however,  if  v  =  0, 


;.<v=o)f  = 


*  * 

s _ 


sinc2^a>Ij 


In  this  case  the  sine  function  becomes 


5inc^<Byj  =  sinc^n5fl)Yj’=S'"c('*^)  =  0  n  =  1,2,3,... 


Thus  all  pixels  are  zero  except  at  co  =  0,  those  exactly  in  the  plane.  For  objects  with  I  vl=  0,  such  as  stars 
in  sidereal  tracking  mode,  the  signal  is  confined  entirely  to  the  to  =  0  plane. 

Expression  ( 1 8)  shows  that  for  v  >  0,  the  signal  falls  off  in  the  to  dimension  with  a  l/n  dependence, 
and  in  fact  a  significant  amount  of  the  signal  power  can  lie  outside  the  primary  signal  plane,  implying 
that  approximations,  assuming  the  signal  is  confined  to  a  single  plane,  can  be  seriously  in  error. 

In  the  spatial  domain  the  energy  per  pixel  [Equation  (4)]  for  I  >p^was 
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(19) 


InBo 


i  Iv 


X 


and  so  the  single  pixel  SNR  is 

2nB  p 

c*  r 


SNR  = 


*  Iv  Ip 

5 

The  matched  filter  operation  integrates  the  entire  energy  of  the  object  for  a  total  signal  strength  of 
2nB  T.  The  total  noise  energy  (above  the  mean)  contained  in  this  streak  is 


so  that  the  output  single  pixel  SNR  from  the  matched  filter  (after  back  transforming)  is 


SNR^ 


=  SNR 

s 


Vq  = 


(20) 


With  this  operation,  one  can  set  a  detection  threshold  a  desired  number  of  standard  deviations  (a) 
above  the  noise  mean.  The  detection  criterion  on  the  filtered  signal  (£^)  is  then  >  a  (E^.  Using  Gaussian 
statistics  the  and  can  now  be  derived.  For  a  Gaussian  distribution  with  standard  deviation  X  = 
the  probability  of  observing  a  noise  contribution  greater  than  aX  is 


oA  a/v2 

where  erfc(x)  is  the  complementary  error  function  normalized  such  that  erfc(O)  =  1 .  The  system  (for 
one  data  set)  is  1  -(1  ~  Th®  probability  of  detecting  a  signal  of  strength 

embedded  in  this  noise  is: 


/’,(£)  =  ! 


E  >aX 

S 


aX-E^ 


E  <ak 

S 
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which  may  be  written 


/■_,(SNR^,ct)=i+i 


£  -aA]  1  1 

V2A  \~  2'^  2 


erf 


SNR^-a 

V2 


SNR^>a  (22a) 


/^(SNR^,  a)  =  — erfc 


aA  -  £  "I  1 
— 7=r^  =~  erfc 
V2A  2 


a-SNRj 

~^r~ 


SNR^  <  a 


(22b) 


Using  Equations  (21)  and  (22),  families  of  receiver  operator  curves  (ROCs)  can  be  generated, 
representing  as  a  function  of  the  filtered  SNR  for  fixed  Pj^.  Because  the  filtered  SNR  is  just  SNRy  = 
SNR^Vft,  one  can  easily  determine  Pj  at  a  given  Pj^  for  a  desired  SNR^  and  velocity,  which  is  shown 
as  a  single  curve  in  Figure  3.  To  use  this,  one  determines  the  input  single  pixel  SNR  and  the  candidate 
velocity,  from  which  one  computes  the  filtered  SNR  (SNRy=  SNR^  -^streak  length ).  Then  for  a  desired 
tlireshold  a  (false  alarm  rate),  the  abscissa  is  determined  and  the  may  then  be  read  off  the  graph. 

In  digital  electronic  implementations,  Fourier  transforms  are  relatively  expensive  operations;  one 
might  ask  if  detection  is  possible  in  the  Fourier  volume,  thereby  eliminating  the  need  for  the  inverse 
transform.  While  the  noise  adds  incoherently,  the  signal  must  add  coherently  with  a  phase  factor  that 
varies  from  point  to  point.  The  calculation  of  this  phase  factor  and  the  subsequent  coherent  addition  is 
nothing  less  than  the  inverse  transform;  thought  of  another  way,  the  noise  density  remains  constant 
through  the  transform;  however,  the  signal  is  transformed  from  a  one-dimensional  line  to  a  two-dimen¬ 
sional  plane.  Thus  the  single  pixel  SNR  is  actually  lower  in  the  Fourier  volume. 


SNR -THRESHOLD 


Figure  3.  General  detection  curve,  fixed  Pj^. 
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3.  RESULTS 


3.1  SIMULATION  SNR  ENHANCEMENT  AND  ¥j^ 

A  three-dimensiona]  matched  filtering  package  was  written  in  C  to  run  under  Unix  and  X- 
Windows/Motif.  The  heart  of  the  package  is  an  n-dimensional  fast  Fourier  transform  algorithm  adapted 
from  Numerical  Recipes  [11].  One  should  be  aware  that  the  results  of  a  fast  fourier  transform  (FFT)  are 
folded  in  frequency  space  as  described  below.  Wells  et  al.  [12]  give  a  good  description  of  the  data 
geometry  involved  in  three-dimensional  transforms.  While  the  ultimate  goal  is  to  process  real  data  sets, 
it  is  important  to  be  able  to  generate  synthetic  data  having  known  properties  to  verify  the  model,  the 
algorithm,  and  its  understanding  under  controlled  circumstances.  Therefore  a  test  pattern-generating  pro¬ 
gram  was  developed  that  generates  three-dimensional  data  sets  composed  of  Gaussian  white  noise  having 
a  specified  mean  and  standard  deviation.  To  this  may  be  added  constant  intensity  objects  moving  at 
specified  velocities  and  randomly  placed  single  pixel  events  corresponding  to  gamma  ray  bombardment 
effects.  The  objects  would  correspond  to  satellites  or  stars  and  have  velocities  in  any  direction  or  zero 
(corresponding  to  stars  in  the  sidereal  case)  with  different  brightness  levels  to  simulate  the  intensity 
difference  between  typical  stars  and  satellites.  The  program  can  generate  image  sets  using  either  a 
5-function  PSF  or  Gaussian  PSF.  Typical  image  sets  using  existing  CCD  cameras  and  electronics  are 
digitized  to  12  bits,  yielding  a  0-  to  4,09S-pixel  gray  scale.  A  typical  noise  mean  is  734  with  a  standard 
deviation  of  7.3.  Stars  can  take  on  values  anywhere  up  to  4,095  and  can,  in  fact,  saturate  the  detector 
causing  several  nonlinear  and  multipixel  effects  that  are  beyond  the  scope  of  this  report.  Similarly, 
satellites  can  have  single  pixel  SNR  values  from  <1  to  >10.  The  single  pixel  SNR  is  measured  according 
to  (value  -  average)/standard  deviation.  The  SNR  of  a  star  is  referred  to  as  the  clutter-to-noise  ratio,  or 
CNR. 

Several  performance  demonstrations  are  shown  using  two-dimensional  images  because  they  are 
easier  to  display  on  paper  and  simple  to  comprehend.  Figure  4  shows  a  two-dimensional  streak  of 
constant  intensity  superimposed  on  a  white  noise  background.  A  matched  filter  was  generated  by  forming 
a  two-dimensional  streak  of  the  same  length  located  at  the  origin  (lower  left)  in  a  field  with  no  noise. 
Figure  5  shows  the  result  of  performing  an  FFT  on  the  image  and  the  filter,  multiplying  the  two,  and  then 
performing  an  inverse  FFT.  The  result  is  the  expected  convolution  of  two  rectangle  functions:  a  triangle 
function.  The  center  of  the  triangle  will  actually  appear  at  the  location  of  one  end  of  the  streak  in  Figure  4, 
but  the  matched  filter  arbitrary  phase  factor  was  used  to  move  the  output  peak  to  the  center  of  the  original 
streak.  If  this  process  is  performed  with  streaks  of  various  lengths,  then  the  data  shown  in  Figure  6,  which 
follows  the  theoretical  curve  quite  well,  is  obtained,  thus  verifying  this  model  and  its  implementation. 

As  the  study  moves  to  full  three-dimensional  data  sets,  one  of  a  number  of  problems  that  arises  is 
how  to  view  and  represent  the  data.  Typically,  each  data  frame  (representing  a  constant  time  slice)  is 
viewed  as  an  x,y  grid.  In  this  case  time  forms  the  third  axis.  One  must  decide  whether  to  represent  these 
data  as  a  series  of  two-dimensional  slices,  as  a  projection,  or  as  a  color  code.  A  maximum  value  pro¬ 
jection  technique  is  well  suited  to  space  surveillance  data  sets  because  it  allows  a  compact  (single  frame) 
representation  without  a  color  requirement.  For  each  {x,y)  point  in  the  three-dimensional  image  volume. 
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Figure  4.  Original  two-dimensional  data:  21  pixel  streak,  single  pixel  SNR  =  8.2. 


Figure  5.  Two-dimensional  matched  filter  output.  Output  SNR  =  36. 
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STREAK  LENGTH  (Pixel) 


Figure  6.  Two-dimensional  matched  fdter  SNR  enhancement.  Input  single  pixel  SNR  =  8.2. 


the  maximum  value  projection  examines  each  data  frame  and  uses  only  the  maximum  value  to  produce 
the  output  two-dimensional  image  frame;  however,  a  strict  maximum  value  projection  technique  suffers 
from  modulation  of  a  constant  intensity  moving  object  because  energy  is  split  between  frames.  For 
example,  an  object  moving  in  the  x  direction  at  1.3  pixels  per  frame  has  a  brightness  in  the  first  pixel 
but  only  1/3  in  the  second.  On  the  next  frame  tte  second  pixel  has  2/3  B^.  Thus  due  to  the  energy 
split  the  second  pixel  shows  an  artificial  reduction,  which  can  be  overcome  by  using  a  pairwise  maximum 
value  projection;  however,  this  technique  is  generally  only  applicable  to  objects  moving  faster  than  1 
pixel  per  frame.  The  effect  of  the  two  techniques  is  shown  in  Figures  7  and  8. 

An  object  translating  parallel  to  the  jc-axis  is  now  used  to  explore  the  Fourier  volume.  Upon  gen¬ 
erating  an  image  volume  and  taking  its  Fourier  transform,  one  finds  (as  predicted)  that  the  Fourier  image 
volume  is  independent  of  the  coordinate.  Taking  a  slice  through  the  k  ,  (O  plane  yields  the  images  shown 
in  Figures  9  and  10.  Two  features  are  worth  noting;  instead  of  a  single  plane  there  are  clearly  two  planes 
containing  signal  in  this  image,  and  the  two  planes  are  not  of  uniform  intensity  but  rather  exhibit  a  roll¬ 
off  phenomenon  along  their  length. 

Figure  9  shows  the  in-plane  roll-off  due  to  a  Gaussian  PSF  having  o  =  0.74  pixels,  which  corre¬ 
sponds  to  one  of  Lincoln's  commonly  used  telescopes.  A  5-function  PSF  transforms  to  yield  a  sine 
function  roll-off  dependent  on  integration  time,  as  shown  in  Figure  10.  The  roll-off  in  this  case  may  be 
written  as  sinc[k  •  v  p,/2],  which  may  be  rewritten  in  one  dimension  as  sinc[j7tv,/A/  ].  The  first  zero  clearly 
occurs  at  pixel  i  =  N/Vj. 
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UNIFORM  INTENSITY  STREAK 


Figure  7.  Sixteen-frame  maximum  value  projection.  Modulation  due  to  energy  split  between  frames. 
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Figure  8.  Sixteen-frame  pairwise  maximum  value  projection. 
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Figure  9.  versus  w  slice  from  Fourier  volume.  Gaussian  P5F  (a  -  0.74  pixels 
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Figure  10.  versus  w  slice  from  Fourier  volume,  Sfunction  PSF. 
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The  signal  decreases  as  sinc[(a)-i-it  •  v^)7/2]as  one  moves  away  from  the  primary  signal  plane. 
Thus  in  a  constant  plane  the  signal  decreases  as  sinc(8/cv^7)/2,  putting  the  first  zero  at  (Sfcv^TV2  =  it.  Using 
dk^  =  IiUNjPj^,  one  can  solve  for  the  pixel  i,  which  contains  the  first  zero: 

p  N  N 

i  =  =  X 

Tv  NVf 

s  t  f 

where  T  =  and  =  P^f^Pt  Yf pixels  per  frame.  The  width  of  the  streak  will  be  2/  pixels. 

The  signal  plane  that  starts  in  the  upper  right  corner  represents  the  negative  frequency  components 
due  to  the  standard  Fourier  space  wraparound  when  the  signal  bandwidth  exceeds  the  sampling  band¬ 
width.  In  the  case  of  a  single  streak,  the  two  signals  remain  separated  at  all  velocities  except  =  1  pixel 
per  frame;  however,  if  multiple  targets  moving  at  different  velocities  are  present,  their  signals  will 
overlap.  The  only  way  to  prevent  this  mixing  of  frequency  components  is  to  assure  that  the  signal  goes 
to  zero  before  reaching  the  Nyquist  frequency  [13]. 

Generating  and  storing  a  large  bank  of  candidate  velocity  filters  covering  the  full  three-dimensional 
image  volume  can  be  prohibitive  for  current  mass  storage  technology  (especially  if  they  are  needed  on¬ 
line  in  a  real-time  processor).  For  example,  a  single  420  x  420  pixel  x  16  frame  image  when  transformed 
to  a  complex  double  precision  (8  bytes)  Fourier  volume  image  requires  45  Mbytes.  Storage  requirements 
can  be  reduced  considerably;  nevertheless  the  data  volume  remains  significant  because  a  typical  filter 
bank  contains  up  to  1,000  filters  [14].  Fortunately,  an  analytical  expression  for  the  Fourier  volume 
representation  of  the  filters  is  available,  thus  enabling  computation  on  the  fly.  The  particularly  simple 
form  wherein  the  signal  is  confined  to  a  single  (albeit  wrapped)  plane  allows  the  Fourier  volume  to  be 
set  equal  to  zero  except  for  the  signal  plane  itself,  which  is  weighted  by  the  sine  function  or  Gaussian 
roll-off  described  above.  The  width  of  this  plane  in  pixels  (derived  above)  is 

Synthetic  three-dimensional  data  sets  were  generated,  containing  noise  and  a  target  moving  at 
various  velocities  and  with  various  constant  intensity  values.  The  complete  matched  filter  operation  was 
applied,  consisting  of  transforming  both  the  input  data  set  and  a  filter  data  set,  multiplying  the  three- 
dimensional  Fourier  volumes,  and  back  transforming  the  product.  As  predicted,  both  the  single  pixel  and 
the  matched  filter  SNR  increased  linearly  with  target  brighmess,  while  the  matched  filter  SNR  increased 
as  the  square  root  of  the  target  streak  length  (for  constant  brightness  per  pixel).  Results  using  the  full 
volume  multiply  and  the  analytically  computed  filter  are  presented  in  Figure  1 1  for  a  constant  brightness 
per  pixel  target  and  in  Figure  12  for  a  constant  brighmess  target.  The  estimates  for  various  thresholds 
were  verified  on  both  original  and  filtered  data.  For  example,  using  64  X  64  pixel  x  16  frames  of  data, 
the  number  of  false  alarms  as  a  function  of  the  threshold  a  was  verified  to  be  n  =  and  X 

could  be  accurately  estimated  from  the  data. 

The  result  of  the  the  full  three-dimensional  matched  filter  is  a  three-dimensional  data  set  in  which 
the  signal  forms  a  triangle  much  as  in  the  two-dimensional  case  (Figure  5),  but  now  it  is  aligned  along 
the  target  trajectory  in  three  dimensions.  A  maximum  value  projection  of  this  signal  would  in  fact 
produce  an  image  similar  to  that  in  Figure  5.  In  comparison,  if  the  predetermined  temporal  plane  con¬ 
taining  the  peak  SNR  were  selected,  then  the  two-dimensional  image  shown  in  Figure  13  would  be 
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Figure  II.  Comparison  of  analytic  matched  filter  with  volume  multiply  for  constant  input  SNR. 
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Figure  12.  SNR  enhancement  as  a  percent  of  theoretical  6-function  PSF. 
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Figure  13.  Two-dimensional  SNR  peak  frame  result. 


obtained.  This  result  is  a  single  frame  of  data  as  opposed  to  the  full  matched  filter  three-dimensional  data 
set,  and  the  signal  peak  is  much  narrower;  it  is  a  triangle  twice  the  width  of  the  single  frame  motion  (2v') 
rather  than  a  triangle  twice  the  width  of  the  entire  data  set  motion  (2vW,).  The  matched  filter  arbitrary 
phase  factor  has  been  used  to  position  the  signal  peak  in  frame  N/2,  frame  8,  for  the  current  16-frame 
data  sets.  This  pitKedure  was  used  to  minimize  wraparound  effects. 

Using  synthetic  images  with  no  jitter  or  Poisson  noise  effects,  the  Fourier  star  slice  removal 
technique  produced  excellent  results.  The  stars  were  completely  eliminated  with  no  residual  artifacts.  The 
impact  on  further  processing,  specifically  the  signal  enhancement  filtering,  was  negligible.  It  was  found 
that  for  all  but  the  lowest  velocity  objects  the  output  SNR  was  not  affected.  At  very  low  velocities  the 
star  and  signal  planes  begin  to  have  a  significant  intersection  volume  leading  to  signal  reduction  in  the 
final  output. 

3,2  REAL  DATA  SNR  ENHANCEMENT  AND  P,  Pj^ 

The  model  assumes  point  objects  moving  with  no  intensity  variation  (signature  or  I^iisson)  and  no 
positional  errors  O't^cr).  In  reality  both  effects  are  present,  meaning  that  one  should  expect  some  discrep¬ 
ancies  between  the  model  and  reality.  For  example,  the  Fourier  domain  star  notch  filter  will  not  com¬ 
pletely  remove  the  stars,  and  the  output  signal  from  a  detected  satellite  will  no  longer  have  the  nice 
triangular  form  seen  in  the  simulated  data.  In  this  section  an  initial  attempt  is  made  to  translate  the 
preceding  theoretical  results  to  real  data. 

An  interesting  and  representative  data  set,  labeled  **geo.dat,”  having  an  average  background  level 
of  709.7  with  a  Gaussian  distribution  noise  of  o  =  9.1,  was  used.  Several  stars  are  present  with  intensities 
ranging  up  to  our  12-bit  limit  of  4,09S.  There  are  two  geosynchronous  satellites  present  in  the  image  set. 


26 


one  relatively  bright  object  and  one  extremely  faint.  For  ease  of  processing  this  data  set  was  clipped  to 
form  two  data  sets  of  64  x  64  pixels  x  16  frames  each.  They  were  labeled  “geob.dat”  to  denote  the  data 
set  containing  the  bright  streak  and  “geow.dat”  to  denote  the  data  set  containing  the  faint  streak.  During 
the  data  collection  process  a  constant  bias  current  is  added  to  the  analog  current  prior  to  the  analog-to-digital 
converter.  To  relate  the  peak  values  of  objects  such  as  stars  to  their  Poisson  noise,  the  bias  must  first  be 
subtracted.  For  the  image  set  geo.dat  this  bias  had  a  digital  value  of  601. 

A  maximum  value  projection  of  the  bright  streak  image  set  is  shown  in  Figure  14,  while  Figure  IS 
is  a  histogram  of  the  noise  distribution.  There  are  five  stars  visible  in  this  image.  Using  a  brighmess  (row, 
column)  notation:  4,095  (56,50);  1,228  (48,51);  2.710  (28,56);  1,164  (22.23);  and  1,124  (37,36).  The 
streak  has  a  maximum  single  pixel  value  of  1,026  and  an  average  value  901 .  Using  the  value  of  901  gives 
an  SNR  =  21.1.  The  streak  profile  also  reveals  variations  (standard  deviation  =  53)  grossly  in  excess  of 
background  noise  level  fluctuations  due  to  frame-to-frame  energy  splitting.  Figure  16  shows  the  Fourier 
volume  (closely  resembling  Figure  9),  the  image  of  which  corresponds  to  a  Gaussian  PSF  (0  =  0.74 
pixels).  In  this  instance  the  star  map  has  been  applied  to  eliminate  the  stars,  followed  by  a  Fourier 
transform. 

An  examination  of  the  stars  reveals  that  the  standard  deviation  computed  for  a  time  sequence  of 
pixels  containing  a  star  is  typically  two  to  three  times  greater  than  predicted  using  -^average  -  bias.  The 
origin  of  these  fluctuations  appears  to  be  atmospheric  or  mount-related  jitter  causing  the  high  gradient 
star  PSFs  to  move  across  pixel  boundaries. 
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Figure  14.  Original  image  data  set,  maximum  value  projection. 
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INTENSITY 


Figure  15.  Nonlinear  least  squares  Gaussian  fit  of  original  data. 


geosb.ftt 


Figure  16.  Fourier  tran^orm  of  geosb.dat  bright  streak.  Star  map  cl  'tter  suppression  applied. 


The  results  of  the  clutter  suppression  algorithms  are  shown  in  Figures  17  through  19.  The  star  slice 
removal  algorithm  reduces  but  does  not  eliminate  the  stars,  while  the  normalizer  produces  a  map  of  single 
pixel  SNR.  CNR  is  computed  using  (4,095-709.7)/  9.1  =  372. 
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geob  rs.inax:  BRIGHT  STREAK 


Figure  17.  Clutter  suppression:  Fourier  star  sticer  applied,  maximum  value  projection. 


g«08b.max:  BRIGHT  STREAK 


Figure  18.  Clutter  suppression:  star  map  mask  applied,  maximum  value  projection. 
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Figure  19.  Clutter  suppression:  maximum  likelihood  normalizer  (F  statistics)  maximum  value  projection. 


Examining  the  original  data  shows  a  streak  length  of  31  pixels,  which  corresponds  to  -1.94 
pixels/frame  and  =  0.  After  applying  the  matched  filter  for  signal  enhancement,  the  final  images  shown 
in  Figures  20  through  22  and  summarized  in  Table  1  are  obtained.  These  figures  represent  the  pairwise 
maximum  value  projection  of  the  full  three-dimensional  result.  In  this  case  the  enhanced  SNR  is  suffi¬ 
ciently  large  that  peak  shifting  due  to  noise  presents  no  problem.  The  two-dimensional  result  obtained 
using  only  the  peak  output  frame  (corresponding  to  the  shift-and-add  result)  is  shown  in  Figure  23  for 
normalizer  clutter  suppression,  corresponding  to  Figure  22.  The  values  for  signal  and  star  peak  represent 
the  image  values  with  the  background  average  subtracted.  The  enhanced  SNR  is  expected  to  be  the  single 
pixel  SNR  multiplied  by  the  square  root  of  the  streak  length:  SNR,^  =  SNR^  ^ 

Both  the  results  in  Table  1  and  Figures  20  through  22  confirm  that  all  three  methods  yield  good  agree¬ 
ment  with  theory;  however,  the  star  slice  clutter  suppression  method  does  allow  some  leakage  (upper  right 
comer  of  Figure  20).  A  binary  threshold  may  be  applied  to  obtain  a  clean  detection  using  any  of  these 
results. 

The  faint  streak  image  set  has  an  average  noise  level  of  713.6  with  a  standard  deviation  of  9.7.  Two 
stars  are  visible  in  the  image:  2108  (20,57)  and  896  (25,60).  Several  significantly  weaker  stars  are  also 
present.  The  streak,  while  present  horizontally  in  the  center  of  the  image,  is  completely  indiscernible.  A 
maximum  value  projection  of  the  original  image  set  is  presented  in  Figure  24.  Figures  25  through  27 
show  the  result  of  the  various  clutter  suppression  techniques  followed  by  the  matched  filter.  In  this  case 
a  series  of  filters  was  used.  The  best  output  SNR  was  obtained  using  v^  =  -1.55  pixels/frame  and  =  0, 
corresponding  to  a  streak  length  of  25  pixels  (see  Table  2). 
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Figure  20.  Fourier  star  slice  +  matched  filter. 


geobjs.  mx2:  BRIGHT  STREAK 


Figure  21.  Star  map  +  Fourier  matched  filter. 
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Figure  22.  Maximum  likelihood  normalizer  -f  Fourier  matched  filter. 


geobn  fs.idt:  BRIGHT  STREAK 


Figure  23.  Peak  output  frame:  maximum  likelihood  normalizer  -f  Fourier  matched  filter. 
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TABLE  1 


geob.dat  Bright  Streak* 


Signal 

Noise 

%  Theoretical 
SNR 

Star 

Peak 

Star 

CNR 

Original  data 

192 

9.1 

21.1 

3386 

372 

Star  map  +  filter 

164 

1.33 

123.2 

106.6 

0 

0 

Star  slice  +  filter 

164 

1.36 

120.4 

103.5 

15.4 

11.3 

Normalizer  -t-  filter 

219 

2.00 

109.7 

94.3 

0 

0 

'Length  =  31  pixels 
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Figure  24.  Original  image  data  set.  maximum  value  projection. 
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geow  rs  fs.mx2:  FAINT  STREAK 


Figure  25.  Fourier  star  slice  +  matched  filter,  pairwise  maximum  value  projection. 


Figure  26.  Star  map  +  Fourier  matched  filter,  pairwise  maximum  value  projection. 
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Figure  27.  Maximum  likelihood  normalizer  +  Fourier  matched  filter,  pairwise  maximum  value  projection. 


TABLE  2 


geow.dat  Paint  Streak* 


Signal 

Noise 

m 

Star 

Peak 

Star 

CNR 

Original  data 

9.66 

1394 

139 

Star  map  +  filter 

10.5 

1.22 

8.5 

0 

0 

Star  slice  +  filter 

10.1 

1.14 

8.8 

12.8 

11.2 

Normalizer  +  filter 

1.39 

0.24 

5.8 

0 

0 

‘Length  =  25  pixels 


As  with  the  bright  streak,  the  normalizer  and  the  star  map  completely  remove  all  the  stars,  including 
many  of  low  intensity,  while  the  Fourier  star  notch  filter  dramatically  reduces  the  star  intensities  but  fails 
to  completely  suppress  them.  In  this  case  the  residual  star  contribution  is  brighter  than  the  enhanced 
streak.  The  original  target  is  too  faint  to  measure  the  SNR  directly;  however,  working  backwards  from 
the  detection  the  unprocessed  SNR  is  estimated  to  be  between  1  and  1.7. 
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To  achieve  detection  a  binary  threshold  (a  =  4)  is  applied  to  the  filter  output  to  obtain  Figures  28 
through  30.  Comparing  the  three  detection  images  it  is  seen  that  using  the  a  priori  threshold,  only  the  star 
map  method  gives  a  clean  detection.  The  star  slice  gives  a  strong  detection  but  also  detects  residual  star 
energy.  The  normalizer  also  shows  a  strong  detection  but  is  accompanied  by  many  random  false  alarms. 
By  experimentally  adjusting  the  threshold  (a  =  3),  the  star  map  and  the  normalizer  give  a  clean  detection 
(Figure  3 1 ).  The  star  slice  technique  has  dramatically  reduced  the  star  contribution,  but  it  remains  the 
brightest  feature  in  the  image  set,  and  thus  a  clean  detection  is  not  possible.  Based  on  Figures  26  and  27, 
while  the  star  map  and  the  normalizer  give  clean  detections,  it  appears  that  the  star  map  gives  a  “stronger” 
detection. 
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Figure  28.  Threshold  image:  Fourier  star  slice  +  matched  filter. 


Figure  29.  Threshold  image:  star  map  +  matched  filler. 
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Figure  30.  Threshold  image:  normalizer  +  matched  filter. 


THRESHOLD  ~  3a 


Figure  31.  Threshold  image:  star  map. 
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4.  SUMMARY 


A  model  of  image  sets  used  in  space  surveillance  applications  has  been  generated  from  which  a 
three-dimensional  Fourier  matched  filter  theory  has  been  developed  and  tested.  SNR  enhancement, 
and  have  been  computed.  Three  clutter  suppression  algorithms  for  star  removal  have  also  been 
examined  and  compared.  Results  using  synthetic  computer-generated  data  and  real  data  agree  extremely 
well  with  theoretical  predictions,  as  do  preliminary  results  using  actual  field  data. 

Although  the  star  slice  method  for  clutter  suppression  does  not  completely  eliminate  the  stars,  it 
can  reduce  them  to  a  level  less  than  all  but  the  faintest  objects;  it  can  also  eliminate  all  objects  moving 
at  a  fixed  velocity,  which  is  an  advantage.  Two  other  star  removal  techniques  completely  eliminate  the 
stars  but  have  added  drawbacks.  A  star  map  mask,  using  a  star  catalog  and  the  telescope  line  of  sight, 
is  computationally  expensive  but  preserves  the  underlying  statistics.  This  approach  achieves  high  perfor¬ 
mance  by  bringing  additional  a  priori  information  to  the  process.  A  maximum  likelihood  normalizer  can 
be  efficiently  implemented,  and  for  stationary  (fixed)  clutter  it  achieves  performance  approaching  that  of 
the  star  map  technique.  An  attractive  feature  of  this  method  is  that  the  normalizer  does  not  use  additional 
information  and  may  thus  be  applied  to  a  wider  range  of  problems. 

Due  to  the  enormous  storage  requirement  associated  with  a  bank  of  precomputed  matched  filters, 
it  is  preferable  (and  in  the  present  application  easy)  to  analytically  compute  the  Fourier  volume  filter  in 
real  time.  By  recognizing  that  the  peak  output  SNR  will  occur  in  a  known  plane,  one  may  reduce  the 
three-dimensional  inverse  Fourier  transform  to  a  two-dimensional  transform,  thus  saving  considerable 
computation. 

The  final  result  is  a  theory  that  quantitatively  describes  many  of  the  features  found  in  typical  space 
surveillance  data  sets.  From  this  theory  and  P^^  predictions  can  be  generated  for  various  cases,  and 
the  relationship  between  camera  noise,  object  brightness,  and  observation  time  required  for  detection  can 
be  estimated.  A  benchmark  of  performance  is  provided  for  comparison  with  various  computationally 
efficient  implementations.  This  work  also  highlights  the  conflict  between  optimizing  detector  configura¬ 
tion  for  maximum  unprocessed  SNR  versus  optimizing  for  the  estimation  of  a  feature.  For  example,  it 
is  desirable  to  have  large  pixels  relative  to  the  PSF  to  maximize  target  signal  in  relation  to  CCD  read¬ 
out  noise.  On  the  other  hand,  the  signal  processing  dictates  pixels  much  smaller  than  the  PSF  to  reduce 
position  dependence  in  the  matched  filter. 
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5.  FUTURE  WORK 


Future  work  should  include  an  analysis  of  a  statistically  significant  number  of  visible-  and  IR-band 
image  sets,  which  would  allow  complete  end-to-end  description  and  verification  of  processing  perfor¬ 
mance  in  terms  of  and  Using  data  sets  collected  in  the  laboratory  with  a  scene  generator  could 
help  clarify  the  distinction  between  Poisson  photon  statistics  and  jitter  effects.  A  discussion  of  the  number 
and  orientation  of  three-dimensional  filters  followed  by  an  automatic  search  routine  would  be  a  required 
step  to  develop  a  deployed  system.  When  using  a  number  of  filters  there  is  some  leakage  of  signal 
between  those  adjacent,  resulting  in  detection  of  a  single  object  at  several  nearby  velocities.  An  analysis 
of  this  phenomenon  and  development  of  suitable  postdetection  algorithms  to  identify  the  single  proper 
velocity  are  needed. 

Following  clutter  suppression  and  signal  enhancement,  the  next  step  in  the  signal  processing  chain 
is  the  clustering  of  detection  points  into  target  streaks  and  the  subsequent  estimation  of  target  location. 
One  unresolved  question  is  the  effect  of  the  signal  enhancement  processing  on  position  estimation,  and 
specifically,  how  will  the  use  of  the  two-dimensional  SNR  peak  rather  than  the  full  three-dimensional 
result  affect  position  estimation? 

Another  useful  direction  would  be  to  determine  the  optimal  approximation  to  matched  filter  pro¬ 
cessing  for  a  specified  level  of  available  computing  power.  Thus  if  matched  filter  processing  is  recog¬ 
nized  as  optimal  but  only  80%  of  the  required  processing  power  is  available,  can  an  optimal  use  of  that 
processing  power  be  determined?  At  least  two  approaches  suggest  themselves.  One  is  to  retain  the 
matched  filter  and  reduce  the  number  of  velocity  filters;  the  other  is  to  design  suboptimal  algorithms  that 
in  some  limit  approach  the  matched  filter. 

Using  the  theory  developed  herein  one  could  discuss  a  number  of  interesting  situations,  one  of 
which  is  the  minimum  resolvable  separation  for  two  objects  on  a  common  trajectory.  A  detector  imple¬ 
mented  on  a  spacecraft  is  subject  to  gamma  ray  bombardment,  leading  to  many  high  intensity,  single  pixel 
events.  Thus  appropriate  gamma  circumvention  techniques  should  be  developed.  Another  area  of  interest 
is  to  estimate  the  minimum  detectable  size  for  low  earth  orbit  debris,  utilizing  the  present  image  acqui¬ 
sition  hardware  and  estimates  of  debris  reflectance.  Finally,  much  of  this  work  should  be  directly  appli¬ 
cable  to  look-down  surveillance  applications  with  radically  different  background  clutter  characteristics. 
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APPENDK 


DISCRETE  FOURIER  TRANSFORM  REPRESENTATION 


A  description  of  the  Fourier  volume  may  be  developed  using  a  discrete  representation  and  a  discrete 
Fourier  transform.  The  physical  intensity  distribution  incident  on  the  focal  plane  is  continuous  in  both 
space  and  time.  A  target  having  a  nonvarying  spatial  distribution  but  allowed  to  vary  in  overall  intensity 
and  translate  with  velocity  v  may  be  written 

I(x,t)  =f(x-x^it))»g(l) 

=  f(x-x^-v^t)»g(t)  .  (A-1) 

The  g(t)  term  allows  intensity  variation  as  a  function  of  time,  such  as  target  signature.  For  this  work 
g(t)  =  1 .  The  starting  position  at  time  r  =  0  is  given  by  x^.  The  function  J(),  representing  the  translating 
target  intensity  distribution,  may  be  represented  by  any  number  of  functions.  For  example. 


/(x -X  -V  t)  = 
'  o  s' 


1 

o 

or 

/  S(x-x  -vt) 

o  o  s 


Next  a  detector  is  assumed  to  have  by  pixels,  each  of  width  X  Wy  It  is  further  assumed 
that  N,  frames  of  data  will  be  collected,  each  having  an  integration  time  T.  The  observed  pixels  will  then 
have  values 


+0* 

f-(*+|)r 

1 - 

-T( 

+ 

1 

H 

1 _ 

8(iJ,k)=  JjJ  f{x-x-^  rect 

2 

rect 

2  •• 

T 

w 

^90 

X 

rect 


y-ii  +  -)w^ 


w 


dxdydt 


(A -2) 


i  =  0,  ...  N  -1 

X 

y  =  0,  ...  N  -1 

J  y 

k  =  0,  ...  N-i 


(x  pixel  index) 
(y  pixel  index) 
(frame  index) 
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which  it  would  be  desirable  to  reduce  to  the  form  =  gO  -  j  -  kv'^- for  reasons 

discussed  below. 

If  the  ^-function  representation  is  used,  then  after  substituting  v  =  v'  w/T  and  z  =  t/T,  the  result  is 


gii,j,k)=  IT  j*rect 


dz 


(A -3) 


The  discrete  Fourier  transform  [15]  is 


N  -1  W  -*  N-\ 


.2nmi  .-2  >017  .2ap4 

G(m,n,p)^-j^  S  X  i  ^  ' 

y  X  y  I  1=0  7=0  k=0 


-  (A-4) 


To  form  a  matched  filter,  the  unknown,  initial  position  information  must  be  separated  from  the  rest 
of  the  expression  [recall  that  in  the  continuous  case.  Equation  (10),  this  process  amounted  to  a  simple 
phase  factor].  Thus  not  only  do  Equations  (A-3)  and  (A-4)  not  lead  to  closed  form  solutions,  they  are 
fundamentally  unsuitable  because  of  this  nonseparable  property.  A  true  matched  filter  has  to  test  not  only 
every  velocity,  but  every  starting  position  for  every  velocity. 

When  instead  of  the  6  function  a  Gaussian  distribution  is  used,  the  result  is  no  better.  In  two 
dimensions  the  problem  is  indeed  tractable,  leading  to 

(*+l)T 

g(i,A:)  =  V2«  al^  I  [erf(((i  +  l)H'-jt^-vr)/<T)-erf((/H'-x^-vf)/cr)]dr 
kT 


which  may  be  evaluated  using 


Jerf(a  +  bjt)<ic 


(a  +  bx)erf(a  +  bx)  + 


I  --ia+bx)^ 

e  ^ 


to  eventually  yield  a  messy  combination  of  error  functions  and  Gaussians  that  are  in  the  form  g{i  -  vk-x^w). 
Unfortunately,  in  three  dimensions  the  problem  has  proven  intractable.  This  formulation  nevertheless  can 
be  used  to  generate  test  data  sets  as  long  as  the  target  velocity  is  confined  parallel  to  one  axis  of  the  pixel 
array.  The  use  of  the  Gaussian  PSF  provides  a  much  closer  representation  of  real  data  sets  than  the  simple 
5-function  form. 

In  light  of  the  foregoing  limitations,  the  model  of  the  satellite  is  used  with  the  following  approxi¬ 
mation: 


(A -5) 
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The  discrete  Fourier  transform  is 
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V  Jr  y  /  /=0  7=0  4=0 

In  this  case  =  InJN^,  dk^  =  liUN^,  d(0  =  2jt/N,  and  to  =  />da).  Several  new  variables 

are  introduced:  z,  =  i-xjw-kv^  and  Zj  =  j-yjw-kv’^,  and  io~ jo~yJ^y'  yielding: 
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(A -6) 


^(2|.Z2^  **  sufficiently  well  behaved  that  periodic  boundary  conditions  can  be  invoked,  the  double 
summation  reduces  to  the  Fourier  transform  of  the  single  frame  target  signature.  Thus, 
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N,  N  N  , 
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where  =  y  = 

This  is  easily  seen  to  be  the  discrete  equivalent  of  Equation  (10).  The  signal  is  maximized  in  the  plane 


p  =  N^ 


mVx  'tVy 
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(A -8) 
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The  matched  filter  is  then  the  complex  conjugate  of  with  =  0: 


H(m,n,p)  = 


Gim,n) 


sin(ff)^^) 


sin(;ry) 


JNNN 
y  jc  y  I 

=  -f— G(m,n)  /cos[;r(A^  -l)y]+isint;r(A('  -I)y]} 

JN_N_.N.  sinOry)  I  '  '  ‘  J 


(A -9) 


y  t 


Unfortunately,  periodic  boundary  conditions  cannot  normally  be  invoked  because  the  signature 
varies  from  frame  to  frame.  It  is  important  to  note  that  this  is  not  an  intrinsic  signature  fIZpZj)  but  one 
that  is  induced  due  to  the  integration  effects  of  the  finite  pixels.  If  the  (Gaussian)  PSF  is  larger  than  a 
single  pixel,  the  signature  will  not  vary  from  frame  to  frame;  however,  in  the  case  where  a  is  on  the  order 
of  a  pixel  size,  significant  velocity-dependent  effects  arise.  The  worst  case  is  a  narrow  PSF  and  low 
velocities,  on  the  order  of  1 .5  pixels  per  frame. 
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SYMBOLS 


i\(t),  Ni(o) 
s(t),  S(co) 
/?(x) 

k.(o 

I(x,t) 

Ws 


Px'Py'Pt 


T1 

n 


Vf,Vf 

G(k,ia) 


probability  of  detection 
probability  of  false  alarm 

number  of  pixels  in  x,y  direction  and  number  of  frames 

Fourier  transform  pair  representing  noise 

Fourier  transform  pair  representing  signal 

autocorrelation  function  (of  noise) 

object  starting  position 

spatial  and  temporal  frequency  components 

focal  plane  intensity  distribution 

star,  noise,  and  signal  contribution  to  focal  plane  intensity 
distribution 

pixel  and  total  data  volume 

pixel  dimensions 

total  data  volume  dimensions 

Gaussian  point  spread  distribution  parameter 

target  velocity 

target  velocity  in  pixels/frame 
target  Gaussian  peak  intensity 
total  signal  energy  per  frame 
signal  power 

average  signal  energy  per  pixel 

observed  noise  standard  deviation  per  pixel 

background  average  value 

total  number  of  pixels  in  a  streak 

noise  energy 

Fourier  volume  per  pixel,  total  Fourier  volume 
Fourier  transform  of  signal 


47 


ij.k 

indices  referring  to  space -tinie  data 

m.n,p 

indices  referring  to  Fourier  volume  data 

a 

threshold  parameter  (in  units  of  standard  deviations) 

SNR,.SNR^ 

signal-to-noise  ratio  of  original  (s)  and  filtered  (/)  data 

X 

PpVQ 

CNR 

clutter-to-noise  ratio  (SNR  of  stars) 

y 

Y=0  defines  Fourier  plane  of  maximum  signal 
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